{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to Brian part 3: Simulations\n",
    "\n",
    "If you haven’t yet read parts 1 and 2 on Neurons and Synapses, go read them first.\n",
    "\n",
    "This tutorial is about managing the slightly more complicated tasks that crop up in research problems, rather than the toy examples we've been looking at so far. So we cover things like inputting sensory data, modelling experimental conditions, etc.\n",
    "\n",
    "As before we start by importing the Brian package and setting up matplotlib for IPython:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from brian2 import *\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multiple runs\n",
    "\n",
    "Let's start by looking at a very common task: doing multiple runs of a simulation with some parameter that changes. Let's start off with something very simple, how does the firing rate of a leaky integrate-and-fire neuron driven by Poisson spiking neurons change depending on its membrane time constant? Let's set that up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# remember, this is here for running separate simulations in the same notebook\n",
    "start_scope() \n",
    "# Parameters\n",
    "num_inputs = 100\n",
    "input_rate = 10*Hz\n",
    "weight = 0.1\n",
    "# Range of time constants\n",
    "tau_range = linspace(1, 10, 30)*ms\n",
    "# Use this list to store output rates\n",
    "output_rates = []\n",
    "# Iterate over range of time constants\n",
    "for tau in tau_range:\n",
    "    # Construct the network each time\n",
    "    P = PoissonGroup(num_inputs, rates=input_rate)\n",
    "    eqs = '''\n",
    "    dv/dt = -v/tau : 1\n",
    "    '''\n",
    "    G = NeuronGroup(1, eqs, threshold='v>1', reset='v=0', method='exact')\n",
    "    S = Synapses(P, G, on_pre='v += weight')\n",
    "    S.connect()\n",
    "    M = SpikeMonitor(G)\n",
    "    # Run it and store the output firing rate in the list\n",
    "    run(1*second)\n",
    "    output_rates.append(M.num_spikes/second)\n",
    "# And plot it\n",
    "plot(tau_range/ms, output_rates)\n",
    "xlabel(r'$\\tau$ (ms)')\n",
    "ylabel('Firing rate (sp/s)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now if you're running the notebook, you'll see that this was a little slow to run. The reason is that for each loop, you're recreating the objects from scratch. We can improve that by setting up the network just once. We store a copy of the state of the network before the loop, and restore it at the beginning of each iteration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_scope() \n",
    "num_inputs = 100\n",
    "input_rate = 10*Hz\n",
    "weight = 0.1\n",
    "tau_range = linspace(1, 10, 30)*ms\n",
    "output_rates = []\n",
    "# Construct the network just once\n",
    "P = PoissonGroup(num_inputs, rates=input_rate)\n",
    "eqs = '''\n",
    "dv/dt = -v/tau : 1\n",
    "'''\n",
    "G = NeuronGroup(1, eqs, threshold='v>1', reset='v=0', method='exact')\n",
    "S = Synapses(P, G, on_pre='v += weight')\n",
    "S.connect()\n",
    "M = SpikeMonitor(G)\n",
    "# Store the current state of the network\n",
    "store()\n",
    "for tau in tau_range:\n",
    "    # Restore the original state of the network\n",
    "    restore()\n",
    "    # Run it with the new value of tau\n",
    "    run(1*second)\n",
    "    output_rates.append(M.num_spikes/second)\n",
    "plot(tau_range/ms, output_rates)\n",
    "xlabel(r'$\\tau$ (ms)')\n",
    "ylabel('Firing rate (sp/s)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's a very simple example of using store and restore, but you can use it in much more complicated situations. For example, you might want to run a long training run, and then run multiple test runs afterwards. Simply put a store after the long training run, and a restore before each testing run.\n",
    "\n",
    "You can also see that the output curve is very noisy and doesn't increase monotonically like we'd expect. The noise is coming from the fact that we run the Poisson group afresh each time. If we only wanted to see the effect of the time constant, we could make sure that the spikes were the same each time (although note that really, you ought to do multiple runs and take an average). We do this by running just the Poisson group once, recording its spikes, and then creating a new `SpikeGeneratorGroup` that will output those recorded spikes each time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_scope() \n",
    "num_inputs = 100\n",
    "input_rate = 10*Hz\n",
    "weight = 0.1\n",
    "tau_range = linspace(1, 10, 30)*ms\n",
    "output_rates = []\n",
    "# Construct the Poisson spikes just once\n",
    "P = PoissonGroup(num_inputs, rates=input_rate)\n",
    "MP = SpikeMonitor(P)\n",
    "# We use a Network object because later on we don't\n",
    "# want to include these objects\n",
    "net = Network(P, MP)\n",
    "net.run(1*second)\n",
    "# And keep a copy of those spikes\n",
    "spikes_i = MP.i\n",
    "spikes_t = MP.t\n",
    "# Now construct the network that we run each time\n",
    "# SpikeGeneratorGroup gets the spikes that we created before\n",
    "SGG = SpikeGeneratorGroup(num_inputs, spikes_i, spikes_t)\n",
    "eqs = '''\n",
    "dv/dt = -v/tau : 1\n",
    "'''\n",
    "G = NeuronGroup(1, eqs, threshold='v>1', reset='v=0', method='exact')\n",
    "S = Synapses(SGG, G, on_pre='v += weight')\n",
    "S.connect()\n",
    "M = SpikeMonitor(G)\n",
    "# Store the current state of the network\n",
    "net = Network(SGG, G, S, M)\n",
    "net.store()\n",
    "for tau in tau_range:\n",
    "    # Restore the original state of the network\n",
    "    net.restore()\n",
    "    # Run it with the new value of tau\n",
    "    net.run(1*second)\n",
    "    output_rates.append(M.num_spikes/second)\n",
    "plot(tau_range/ms, output_rates)\n",
    "xlabel(r'$\\tau$ (ms)')\n",
    "ylabel('Firing rate (sp/s)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see that now there is much less noise and it increases monotonically because the input spikes are the same each time, meaning we're seeing the effect of the time constant, not the random spikes.\n",
    "\n",
    "Note that in the code above, we created `Network` objects. The reason is that in the loop, if we just called `run` it would try to simulate all the objects, including the Poisson neurons ``P``, and we only want to run that once. We use `Network` to specify explicitly which objects we want to include.\n",
    "\n",
    "The techniques we've looked at so far are the conceptually most simple way to do multiple runs, but not always the most efficient. Since there's only a single output neuron in the model above, we can simply duplicate that output neuron and make the time constant a parameter of the group."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_scope() \n",
    "num_inputs = 100\n",
    "input_rate = 10*Hz\n",
    "weight = 0.1\n",
    "tau_range = linspace(1, 10, 30)*ms\n",
    "num_tau = len(tau_range)\n",
    "P = PoissonGroup(num_inputs, rates=input_rate)\n",
    "# We make tau a parameter of the group\n",
    "eqs = '''\n",
    "dv/dt = -v/tau : 1\n",
    "tau : second\n",
    "'''\n",
    "# And we have num_tau output neurons, each with a different tau\n",
    "G = NeuronGroup(num_tau, eqs, threshold='v>1', reset='v=0', method='exact')\n",
    "G.tau = tau_range\n",
    "S = Synapses(P, G, on_pre='v += weight')\n",
    "S.connect()\n",
    "M = SpikeMonitor(G)\n",
    "# Now we can just run once with no loop\n",
    "run(1*second)\n",
    "output_rates = M.count/second # firing rate is count/duration\n",
    "plot(tau_range/ms, output_rates)\n",
    "xlabel(r'$\\tau$ (ms)')\n",
    "ylabel('Firing rate (sp/s)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see that this is much faster again! It's a little bit more complicated conceptually, and it's not always possible to do this trick, but it can be much more efficient if it's possible.\n",
    "\n",
    "Let's finish with this example by having a quick look at how the mean and standard deviation of the interspike intervals depends on the time constant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trains = M.spike_trains()\n",
    "isi_mu = full(num_tau, nan)*second\n",
    "isi_std = full(num_tau, nan)*second\n",
    "for idx in range(num_tau):\n",
    "    train = diff(trains[idx])\n",
    "    if len(train)>1:\n",
    "        isi_mu[idx] = mean(train)\n",
    "        isi_std[idx] = std(train)\n",
    "errorbar(tau_range/ms, isi_mu/ms, yerr=isi_std/ms)\n",
    "xlabel(r'$\\tau$ (ms)')\n",
    "ylabel('Interspike interval (ms)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that we used the ``spike_trains()`` method of `SpikeMonitor`. This is a dictionary with keys being the indices of the neurons and values being the array of spike times for that neuron.\n",
    "\n",
    "## Changing things during a run\n",
    "\n",
    "Imagine an experiment where you inject current into a neuron, and change the amplitude randomly every 10 ms. Let's see if we can model that using a Hodgkin-Huxley type neuron."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_scope()\n",
    "# Parameters\n",
    "area = 20000*umetre**2\n",
    "Cm = 1*ufarad*cm**-2 * area\n",
    "gl = 5e-5*siemens*cm**-2 * area\n",
    "El = -65*mV\n",
    "EK = -90*mV\n",
    "ENa = 50*mV\n",
    "g_na = 100*msiemens*cm**-2 * area\n",
    "g_kd = 30*msiemens*cm**-2 * area\n",
    "VT = -63*mV\n",
    "# The model\n",
    "eqs_HH = '''\n",
    "dv/dt = (gl*(El-v) - g_na*(m*m*m)*h*(v-ENa) - g_kd*(n*n*n*n)*(v-EK) + I)/Cm : volt\n",
    "dm/dt = 0.32*(mV**-1)*(13.*mV-v+VT)/\n",
    "    (exp((13.*mV-v+VT)/(4.*mV))-1.)/ms*(1-m)-0.28*(mV**-1)*(v-VT-40.*mV)/\n",
    "    (exp((v-VT-40.*mV)/(5.*mV))-1.)/ms*m : 1\n",
    "dn/dt = 0.032*(mV**-1)*(15.*mV-v+VT)/\n",
    "    (exp((15.*mV-v+VT)/(5.*mV))-1.)/ms*(1.-n)-.5*exp((10.*mV-v+VT)/(40.*mV))/ms*n : 1\n",
    "dh/dt = 0.128*exp((17.*mV-v+VT)/(18.*mV))/ms*(1.-h)-4./(1+exp((40.*mV-v+VT)/(5.*mV)))/ms*h : 1\n",
    "I : amp\n",
    "'''\n",
    "group = NeuronGroup(1, eqs_HH,\n",
    "                    threshold='v > -40*mV',\n",
    "                    refractory='v > -40*mV',\n",
    "                    method='exponential_euler')\n",
    "group.v = El\n",
    "statemon = StateMonitor(group, 'v', record=True)\n",
    "spikemon = SpikeMonitor(group, variables='v')\n",
    "figure(figsize=(9, 4))\n",
    "for l in range(5):\n",
    "    group.I = rand()*50*nA\n",
    "    run(10*ms)\n",
    "    axvline(l*10, ls='--', c='k')\n",
    "axhline(El/mV, ls='-', c='lightgray', lw=3)\n",
    "plot(statemon.t/ms, statemon.v[0]/mV, '-b')\n",
    "plot(spikemon.t/ms, spikemon.v/mV, 'ob')\n",
    "xlabel('Time (ms)')\n",
    "ylabel('v (mV)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the code above, we used a loop over multiple runs to achieve this. That's fine, but it's not the most efficient way to do it because each time we call ``run`` we have to do a lot of initialisation work that slows everything down. It also won't work as well with the more efficient standalone mode of Brian. Here's another way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_scope()\n",
    "group = NeuronGroup(1, eqs_HH,\n",
    "                    threshold='v > -40*mV',\n",
    "                    refractory='v > -40*mV',\n",
    "                    method='exponential_euler')\n",
    "group.v = El\n",
    "statemon = StateMonitor(group, 'v', record=True)\n",
    "spikemon = SpikeMonitor(group, variables='v')\n",
    "# we replace the loop with a run_regularly\n",
    "group.run_regularly('I = rand()*50*nA', dt=10*ms)\n",
    "run(50*ms)\n",
    "figure(figsize=(9, 4))\n",
    "# we keep the loop just to draw the vertical lines\n",
    "for l in range(5):\n",
    "    axvline(l*10, ls='--', c='k')\n",
    "axhline(El/mV, ls='-', c='lightgray', lw=3)\n",
    "plot(statemon.t/ms, statemon.v[0]/mV, '-b')\n",
    "plot(spikemon.t/ms, spikemon.v/mV, 'ob')\n",
    "xlabel('Time (ms)')\n",
    "ylabel('v (mV)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We've replaced the loop that had multiple ``run`` calls with a ``run_regularly``. This makes the specified block of code run every ``dt=10*ms``. The ``run_regularly`` lets you run code specific to a single `NeuronGroup`, but sometimes you might need more flexibility. For this, you can use `network_operation` which lets you run arbitrary Python code (but won't work with the standalone mode)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_scope()\n",
    "group = NeuronGroup(1, eqs_HH,\n",
    "                    threshold='v > -40*mV',\n",
    "                    refractory='v > -40*mV',\n",
    "                    method='exponential_euler')\n",
    "group.v = El\n",
    "statemon = StateMonitor(group, 'v', record=True)\n",
    "spikemon = SpikeMonitor(group, variables='v')\n",
    "# we replace the loop with a network_operation\n",
    "@network_operation(dt=10*ms)\n",
    "def change_I():\n",
    "    group.I = rand()*50*nA\n",
    "run(50*ms)\n",
    "figure(figsize=(9, 4))\n",
    "for l in range(5):\n",
    "    axvline(l*10, ls='--', c='k')\n",
    "axhline(El/mV, ls='-', c='lightgray', lw=3)\n",
    "plot(statemon.t/ms, statemon.v[0]/mV, '-b')\n",
    "plot(spikemon.t/ms, spikemon.v/mV, 'ob')\n",
    "xlabel('Time (ms)')\n",
    "ylabel('v (mV)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's extend this example to run on multiple neurons, each with a different capacitance to see how that affects the behaviour of the cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_scope()\n",
    "N = 3\n",
    "eqs_HH_2 = '''\n",
    "dv/dt = (gl*(El-v) - g_na*(m*m*m)*h*(v-ENa) - g_kd*(n*n*n*n)*(v-EK) + I)/C : volt\n",
    "dm/dt = 0.32*(mV**-1)*(13.*mV-v+VT)/\n",
    "    (exp((13.*mV-v+VT)/(4.*mV))-1.)/ms*(1-m)-0.28*(mV**-1)*(v-VT-40.*mV)/\n",
    "    (exp((v-VT-40.*mV)/(5.*mV))-1.)/ms*m : 1\n",
    "dn/dt = 0.032*(mV**-1)*(15.*mV-v+VT)/\n",
    "    (exp((15.*mV-v+VT)/(5.*mV))-1.)/ms*(1.-n)-.5*exp((10.*mV-v+VT)/(40.*mV))/ms*n : 1\n",
    "dh/dt = 0.128*exp((17.*mV-v+VT)/(18.*mV))/ms*(1.-h)-4./(1+exp((40.*mV-v+VT)/(5.*mV)))/ms*h : 1\n",
    "I : amp\n",
    "C : farad\n",
    "'''\n",
    "group = NeuronGroup(N, eqs_HH_2,\n",
    "                    threshold='v > -40*mV',\n",
    "                    refractory='v > -40*mV',\n",
    "                    method='exponential_euler')\n",
    "group.v = El\n",
    "# initialise with some different capacitances\n",
    "group.C = array([0.8, 1, 1.2])*ufarad*cm**-2*area\n",
    "statemon = StateMonitor(group, variables=True, record=True)\n",
    "# we go back to run_regularly\n",
    "group.run_regularly('I = rand()*50*nA', dt=10*ms)\n",
    "run(50*ms)\n",
    "figure(figsize=(9, 4))\n",
    "for l in range(5):\n",
    "    axvline(l*10, ls='--', c='k')\n",
    "axhline(El/mV, ls='-', c='lightgray', lw=3)\n",
    "plot(statemon.t/ms, statemon.v.T/mV, '-')\n",
    "xlabel('Time (ms)')\n",
    "ylabel('v (mV)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So that runs, but something looks wrong! The injected currents look like they're different for all the different neurons! Let's check:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(statemon.t/ms, statemon.I.T/nA, '-')\n",
    "xlabel('Time (ms)')\n",
    "ylabel('I (nA)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sure enough, it's different each time. But why? We wrote ``group.run_regularly('I = rand()*50*nA', dt=10*ms)`` which seems like it should give the same value of I for each neuron. But, like threshold and reset statements, ``run_regularly`` code is interpreted as being run separately for each neuron, and because I is a parameter, it can be different for each neuron. We can fix this by making I into a *shared* variable, meaning it has the same value for each neuron."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_scope()\n",
    "N = 3\n",
    "eqs_HH_3 = '''\n",
    "dv/dt = (gl*(El-v) - g_na*(m*m*m)*h*(v-ENa) - g_kd*(n*n*n*n)*(v-EK) + I)/C : volt\n",
    "dm/dt = 0.32*(mV**-1)*(13.*mV-v+VT)/\n",
    "    (exp((13.*mV-v+VT)/(4.*mV))-1.)/ms*(1-m)-0.28*(mV**-1)*(v-VT-40.*mV)/\n",
    "    (exp((v-VT-40.*mV)/(5.*mV))-1.)/ms*m : 1\n",
    "dn/dt = 0.032*(mV**-1)*(15.*mV-v+VT)/\n",
    "    (exp((15.*mV-v+VT)/(5.*mV))-1.)/ms*(1.-n)-.5*exp((10.*mV-v+VT)/(40.*mV))/ms*n : 1\n",
    "dh/dt = 0.128*exp((17.*mV-v+VT)/(18.*mV))/ms*(1.-h)-4./(1+exp((40.*mV-v+VT)/(5.*mV)))/ms*h : 1\n",
    "I : amp (shared) # everything is the same except we've added this shared\n",
    "C : farad\n",
    "'''\n",
    "group = NeuronGroup(N, eqs_HH_3,\n",
    "                    threshold='v > -40*mV',\n",
    "                    refractory='v > -40*mV',\n",
    "                    method='exponential_euler')\n",
    "group.v = El\n",
    "group.C = array([0.8, 1, 1.2])*ufarad*cm**-2*area\n",
    "statemon = StateMonitor(group, 'v', record=True)\n",
    "group.run_regularly('I = rand()*50*nA', dt=10*ms)\n",
    "run(50*ms)\n",
    "figure(figsize=(9, 4))\n",
    "for l in range(5):\n",
    "    axvline(l*10, ls='--', c='k')\n",
    "axhline(El/mV, ls='-', c='lightgray', lw=3)\n",
    "plot(statemon.t/ms, statemon.v.T/mV, '-')\n",
    "xlabel('Time (ms)')\n",
    "ylabel('v (mV)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ahh, that's more like it!\n",
    "\n",
    "## Adding input\n",
    "\n",
    "Now let's think about a neuron being driven by a sinusoidal input. Let's go back to a leaky integrate-and-fire to simplify the equations a bit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_scope()\n",
    "A = 2.5\n",
    "f = 10*Hz\n",
    "tau = 5*ms\n",
    "eqs = '''\n",
    "dv/dt = (I-v)/tau : 1\n",
    "I = A*sin(2*pi*f*t) : 1\n",
    "'''\n",
    "G = NeuronGroup(1, eqs, threshold='v>1', reset='v=0', method='euler')\n",
    "M = StateMonitor(G, variables=True, record=True)\n",
    "run(200*ms)\n",
    "plot(M.t/ms, M.v[0], label='v')\n",
    "plot(M.t/ms, M.I[0], label='I')\n",
    "xlabel('Time (ms)')\n",
    "ylabel('v')\n",
    "legend(loc='best');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far, so good and the sort of thing we saw in the first tutorial. Now, what if that input current were something we had recorded and saved in a file? In that case, we can use `TimedArray`. Let's start by reproducing the picture above but using `TimedArray`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_scope()\n",
    "A = 2.5\n",
    "f = 10*Hz\n",
    "tau = 5*ms\n",
    "# Create a TimedArray and set the equations to use it\n",
    "t_recorded = arange(int(200*ms/defaultclock.dt))*defaultclock.dt\n",
    "I_recorded = TimedArray(A*sin(2*pi*f*t_recorded), dt=defaultclock.dt)\n",
    "eqs = '''\n",
    "dv/dt = (I-v)/tau : 1\n",
    "I = I_recorded(t) : 1\n",
    "'''\n",
    "G = NeuronGroup(1, eqs, threshold='v>1', reset='v=0', method='exact')\n",
    "M = StateMonitor(G, variables=True, record=True)\n",
    "run(200*ms)\n",
    "plot(M.t/ms, M.v[0], label='v')\n",
    "plot(M.t/ms, M.I[0], label='I')\n",
    "xlabel('Time (ms)')\n",
    "ylabel('v')\n",
    "legend(loc='best');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that for the example where we put the ``sin`` function directly in the equations, we had to use the ``method='euler'`` argument because the exact integrator wouldn't work here (try it!). However, ``TimedArray`` is considered to be constant over its time step and so the linear integrator can be used. This means you won't get the same behaviour from these two methods for two reasons. Firstly, the numerical integration methods ``exact`` and ``euler`` give slightly different results. Secondly, ``sin`` is not constant over a timestep whereas ``TimedArray`` is.\n",
    "\n",
    "Now just to show that ``TimedArray`` works for arbitrary currents, let's make a weird \"recorded\" current and run it on that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_scope()\n",
    "A = 2.5\n",
    "f = 10*Hz\n",
    "tau = 5*ms\n",
    "# Let's create an array that couldn't be\n",
    "# reproduced with a formula\n",
    "num_samples = int(200*ms/defaultclock.dt)\n",
    "I_arr = zeros(num_samples)\n",
    "for _ in range(100):\n",
    "    a = randint(num_samples)\n",
    "    I_arr[a:a+100] = rand()\n",
    "I_recorded = TimedArray(A*I_arr, dt=defaultclock.dt)\n",
    "eqs = '''\n",
    "dv/dt = (I-v)/tau : 1\n",
    "I = I_recorded(t) : 1\n",
    "'''\n",
    "G = NeuronGroup(1, eqs, threshold='v>1', reset='v=0', method='exact')\n",
    "M = StateMonitor(G, variables=True, record=True)\n",
    "run(200*ms)\n",
    "plot(M.t/ms, M.v[0], label='v')\n",
    "plot(M.t/ms, M.I[0], label='I')\n",
    "xlabel('Time (ms)')\n",
    "ylabel('v')\n",
    "legend(loc='best');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let's finish on an example that actually reads in some data from a file. See if you can work out how this example works."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_scope()\n",
    "from matplotlib.image import imread\n",
    "img = (1-imread('brian.png'))[::-1, :, 0].T\n",
    "num_samples, N = img.shape\n",
    "ta = TimedArray(img, dt=1*ms) # 228\n",
    "A = 1.5\n",
    "tau = 2*ms\n",
    "eqs = '''\n",
    "dv/dt = (A*ta(t, i)-v)/tau+0.8*xi*tau**-0.5 : 1\n",
    "'''\n",
    "G = NeuronGroup(N, eqs, threshold='v>1', reset='v=0', method='euler')\n",
    "M = SpikeMonitor(G)\n",
    "run(num_samples*ms)\n",
    "plot(M.t/ms, M.i, '.k', ms=3)\n",
    "xlim(0, num_samples)\n",
    "ylim(0, N)\n",
    "xlabel('Time (ms)')\n",
    "ylabel('Neuron index');"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
